C
      SUBROUTINE GRIDMP (MAPSET)
C
C     GRIDMP COMPUTES THE INTERPOLATION GRID INTERSECTIONS AND THE
C     INTERPOLATION ERRORS IN IMAGE PIXEL LOCATIONS
C     ALLOWANCE IS MADE FOR THE STARTING NORTHING VALUE GREATER THAN
C     THE STOPPING VALUE (AS IN THE UTM SYSTEM)
C     GRID IS COORDINATES (UTM & PIXEL) OF HORIZONTAL & VERTICAL GRID
C     INTERSECTIONS
C
      DIMENSION GRIDEN(2,30), GRIDSL(2,30,30), NUMPXC(30), CPXVAL(30),
     .NUMLNC(30), ERRCEL(29)
      LOGICAL*1 CELLCG(29,29)
      REAL*4 MVPOFF, MVPINV, MVP, NORTH, NCEN
      LOGICAL MAPSET
      COMMON /PIXDAT/ NLINPX, NPIXLN, NHGL, NVGL, SRATE, NLPI, NSPL,
     .                KDUM, INTERP, NB
      COMMON /MAPDAT/ STRTNO, STOPNO, STRTEA, STOPEA
      COMMON /  GRID/ GRIDEN, GRIDSL, NUMPXC, CPXVAL, NUMLNC, CELLCG
C
C     COMPUTE UTM COORDINATES AT IMAGE CORNERS
      IF (MAPSET) GO TO 900
      FLINPX=NLINPX
      FPIXLN=NPIXLN
      MVP = MVPOFF(1.0)
      CRV = ERCURV(1.0-MVP)
      PS = 1.0 - DELAY(1) - MVP - CRV
      CALL EVPOLY (3, PS, 1.0, AE)
      CALL EVPOLY (4, PS, 1.0, AN)
      PS = 1.0 - DELAY(NLINPX) - MVP - CRV
      CALL EVPOLY (3, PS, FLINPX, CE)
      CALL EVPOLY (4, PS, FLINPX, CN)
      MVP = MVPOFF(FPIXLN)
      CRV = ERCURV(FPIXLN-MVP)
      PS = FPIXLN - DELAY(1) - MVP - CRV
      CALL EVPOLY (3, PS, 1.0, BE)
      CALL EVPOLY (4, PS, 1.0, BN)
      PS = FPIXLN - DELAY(NLINPX) - MVP - CRV
      CALL EVPOLY (3, PS, FLINPX, DE)
      CALL EVPOLY (4, PS, FLINPX, DN)
      IF (AN.LT.CN) STRTNO = AMIN1 (AN, BN)
      IF (AN.GT.CN) STRTNO = AMAX1 (AN, BN)
      IF (AN.LT.CN) STOPNO = AMAX1 (CN, DN)
      IF (AN.GT.CN) STOPNO = AMIN1 (CN, DN)
      STRTEA = AMIN1 (AE, CE)
      STOPEA = AMAX1 (BE, DE)
C
C     COMPUTE GX AND GY GRID SPACING AND OUTPUT IMAGE SIZE
  900 NHGL1=NHGL-1
      NVGL1=NVGL-1
      GX = (STOPEA-STRTEA) / NVGL1
      GY = (STOPNO-STRTNO) / NHGL1
      NSPL = (STOPEA-STRTEA+SRATE) / SRATE + 0.5
      NLPI = (ABS(STOPNO-STRTNO)+SRATE) / SRATE + 0.5
C
C     CALCULATE GRID INTERSECTIONS BY STEPPING THROUGH UTM GRID AND
C     TRANSFORMING TO IMAGE (PIXEL) COORDINATES
        DO 1010 I=1,NVGL
        EAST = STRTEA + (I-1)*GX
        GRIDEN(1,I) = EAST
          DO 1000 J=1,NHGL
          NORTH = STRTNO + (J-1)*GY
          GRIDEN(2,J) = NORTH
          CALL EVPOLY (1, EAST, NORTH, PS)
          CALL EVPOLY (2, EAST, NORTH, PL)
C
C         ADD INVERSES OF MIRROR VELOCITY PROFILE AND EARTH CURVATURE
C         CORRECTIONS TO GET RAW PIXEL COORDINATES AT GRID POINTS
          LINE = PL + SIGN(0.5,PL)
          PSCCT = PS + DELAY(LINE)
          MVP = MVPINV(PSCCT)
          CRV = ERCURI(PSCCT+MVP)
          PS = PS + MVP + CRV
          GRIDSL(1,J,I) = PS
          GRIDSL(2,J,I) = PL
 1000     CONTINUE
 1010   CONTINUE
      WRITE (6,1050)
      DO 1011 J=1,NHGL
 1011 WRITE (6,1051) (GRIDSL(1,J,K), GRIDSL(2,J,K), GRIDEN(1,K),
     .GRIDEN(2,J), K=1,NVGL)
C
C       COMPUTE FINAL GRID ERRORS
          WRITE (6,1065)
        ERRMAX=0.
        DO 1030 I4=2,NHGL
        I41=I4-1
          DO 1020 J4=2,NVGL
          J41=J4-1
C
C         COMPUTE PIXEL COORDINATES BY BILINEAR INTERPOLATION IN A CELL
          PSB = (GRIDSL(1,I4,J4)+GRIDSL(1,I41,J4)+GRIDSL(1,I41,J41)
     .    +GRIDSL(1,I4,J41)) / 4.0
          PLB = (GRIDSL(2,I4,J4)+GRIDSL(2,I41,J4)+GRIDSL(2,I41,J41)
     .    +GRIDSL(2,I4,J41)) / 4.0
C
C         COMPUTE CENTER EASTING AND NORTHING
          ECEN = (GRIDEN(1,J4)+GRIDEN(1,J41)) / 2.0
          NCEN = (GRIDEN(2,I4)+GRIDEN(2,I41)) / 2.0
C
C         COMPUTE PIXEL COORDINATES BY FUNCTION
          CALL EVPOLY (1, ECEN, NCEN, PS)
          CALL EVPOLY (2, ECEN, NCEN, PL)
C
C         ADD INVERSES OF MIRROR VELOCITY PROFILE AND EARTH CURVATURE
C         CORRECTIONS TO GRID POINTS
          LINE = PL + SIGN(0.5,PL)
          PSCCT = PS + DELAY(LINE)
          MVP = MVPINV(PSCCT)
          CRV = ERCURI(PSCCT+MVP)
          PS = PS + MVP + CRV
C
C         COMPUTE ERROR
          ERR=SQRT((PSB-PS)**2+(PLB-PL)**2)
          ERRCEL(J41) = ERR
          ERRMAX = AMAX1 (ERRMAX, ERR)
 1020     CONTINUE
        WRITE (6,1070) (ERRCEL(J), J=1,NVGL1)
 1030   CONTINUE
      WRITE (6,1060) ERRMAX
C
C     COMPUTE TABLES OF FIRST OUTPUT SAMPLE NUMBERS AND COORDINATES AND
C     OUTPUT LINE NUMBERS IN EACH CELL
      NUMPXC(1) = 1
      CPXVAL(1) = GRIDEN(1,1)
      NUMLNC(1) = 1
      DO 1038 J=2,NVGL1
      NUMPXC(J) = (GRIDEN(1,J)-GRIDEN(1,1)+SRATE)/SRATE + 1.0
 1038 CPXVAL(J) = GRIDEN(1,1) + (NUMPXC(J)-1)*SRATE
      DO 1039 J=2,NHGL1
 1039 NUMLNC(J) = (ABS(GRIDEN(2,J)-GRIDEN(2,1))+SRATE)/SRATE + 1.0
      NUMPXC(NVGL) = NSPL + 1
      NUMLNC(NHGL) = NLPI + 1
C
      WRITE (6,1052)
      DO 1045 J=1,29
      IF (J.LE.NHGL1.OR.J.LE.NVGL1) WRITE (6,1054) J
      IF (J.LE.NVGL1) WRITE (6,1055) NUMPXC(J), CPXVAL(J)
 1045 IF (J.LE.NHGL1) WRITE (6,1056) NUMLNC(J)
      WRITE (6,1053) SRATE, NLPI, NSPL
      RETURN
C
 1050 FORMAT ('1',30X,'INTERPOLATION GRID INTERSECTIONS'/31X,32('*')//
     .14X,'SAMPLE',16X,'LINE',13X,'EASTING',12X,'NORTHING')
 1051 FORMAT (/(4F20.3))
 1052 FORMAT ('1'/30X,'GRID CELL STARTING POINTS'/30X,25('*')//10X,'CELL
     . NO.',5X,'SAMPLE NO.',5X,'COORD. VALUE',5X,'LINE NO.'/)
 1053 FORMAT (///21X,'OUTPUT PIXEL SPACING =',F7.3,'  UNITS'/21X,'OUTPUT
     . IMAGE SIZE =',I5,' LINES BY',I5,' SAMPLES')
 1054 FORMAT (I15)
 1055 FORMAT ('+',I30,F18.3)
 1056 FORMAT ('+',I60)
 1060 FORMAT (//21X,'MAXIMUM GRID ERROR =',1PE10.2,'  PIXELS')
 1065 FORMAT ('1',20X,'MAGNITUDE OF INTERPOLATION ERROR AT GRID CELL CEN
     .TERS'/21X,53('*'))
 1070 FORMAT (/(1P14E9.2))
      END
C
      SUBROUTINE CELLMP (CORE)
C
C     *COMPUTE DIMENSIONS OF INPUT ARRAY WHICH CAN BE HELD IN CORE AT A
C      TIME AND NUMBER OF PARTITIONS REQUIRED.
C     PIX AND ROW VARIABLES HAVE DIRECTIONS APPENDED, E.G. PIXNW IS
C     NORTHWEST CORNER PIXEL OF A CELL
C
      INTEGER ROWNW, PIXNW, ROWNE, PIXNE, ROWSE, PIXSE, ROWSW, PIXSW,
     .RN, RS, SW, SE, R1, R2, S1, S2, CORE, RSTRT
      DIMENSION GRIDEN(2,30), GRIDSL(2,30,30), NUMPXC(30), CPXVAL(30),
     .NUMLNC(30)
      LOGICAL*1 CELLCG(29,29)
      COMMON /PIXDAT/ NLINPX, NPIXLN, NHGL, NVGL, SRATE, NLPI, NSPL,
     .                KDUM, INTRP, NB
      COMMON /  GRID/ GRIDEN, GRIDSL, NUMPXC, CPXVAL, NUMLNC, CELLCG
      COMMON /SEGMNT/ NCGP, NBSAM(10), NBLIN(10), NSTRT(10),
     .                RSTRT(29,10)
C
C     INITIALIZE CELL COLUMN GROUP INDICATOR (CELLCG)
      NHGL1 = NHGL - 1
      NVGL1 = NVGL - 1
      DO 1100 I=1,NHGL1
 1100 CELLCG(I,1) = 222
C
C     FIND NCGP, THE NUMBER OF COLUMN GROUPS.
C     S1 IS FIRST INPUT PIXEL FOR THIS COLUMN GROUP
      NCGP = 0
      S1 = 1
C
C     START A NEW COLUMN GROUP OR SEGMENT
 1800 CONTINUE
      NCGP = NCGP + 1
      MINNC = 1000000
C
C     LOOP OVER CELL ROWS
      DO 2010 I=1,NHGL1
      R1 = 1000000
      R2 = 0
      S2 = 0
C
C       FIND FIRST CELL IN THE ROW AVAILABLE FOR THIS COLUMN GROUP AND
C       START THE LOOP THERE (CELL IC1)
        DO 2001 IC1=1,NVGL1
        IF (CELLCG(I,IC1).EQ.222) GO TO 2002
 2001   CONTINUE
        GO TO 2010
C
C         FIND RANGE OF IMAGE COORDINATES REQUIRED FOR THE OUTPUT CELL
 2002     DO 2000 J=IC1,NVGL1
          ROWNW = GRIDSL(2,I,J) - 1.0
          PIXNW = GRIDSL(1,I,J) + DELAY(ROWNW) - 2.0
          ROWNE = GRIDSL(2,I,J+1) - 1.0
          PIXNE = GRIDSL(1,I,J+1) + DELAY(ROWNE) + 3.0
          ROWSE = GRIDSL(2,I+1,J+1) + 2.0
          PIXSE = GRIDSL(1,I+1,J+1) + DELAY(ROWSE) + 3.0
          ROWSW = GRIDSL(2,I+1,J) + 2.0
          PIXSW = GRIDSL(1,I+1,J) + DELAY(ROWSW) - 2.0
          RN = MAX0 (MIN0(ROWNW,ROWNE), 1)
          RS = MIN0 (MAX0(ROWSW,ROWSE), NLINPX)
          SW = MAX0 (MIN0(PIXNW,PIXSW), 1)
          SE = MIN0 (MAX0(PIXNE,PIXSE), NPIXLN)
C
C         CHECK FOR IMAGE COORDINATES OUTSIDE THE INPUT IMAGE
          IF (RN.GT.NLINPX.OR.RS.LT.1.OR.SW.GT.NPIXLN.OR.SE.LT.1)
     .    GO TO 1920
C
C           FIND RANGE OF IMAGE COORDINATES FOR THIS ROW OF CELLS AND
C           FIND NO. OF OUTPUT CELLS THAT CAN BE FILLED GIVEN CORE LIMIT
            R1 = MIN0 (R1, RN)
            R2 = MAX0 (R2, RS)
            S2 = MAX0 (S2, SE)
            NR = R2 - R1 + 1
            NC = S2 - S1 + 1
            IF (NR*NC.LE.CORE) GO TO 2005
C
C           CORE REQUIRED EXCEEDS THAT AVAILABLE
C           GO TO NEXT ROW OF CELLS
            MINNC = MIN0 (NCSAVE, MINNC)
            CELLCG(I,J) = 255
            GO TO 2010
C
 1920     CELLCG(I,J) = 0
          GO TO 2000
C
C           CORE IS AVAILABLE.  SAVE REQUIRED DIMENSIONS.
 2005       NCSAVE = NC
            CELLCG(I,J) = 222
 2000       CONTINUE
 2010 CONTINUE
C
      LASTPX = S1 + MINNC - 1
      IF (LASTPX.LT.NPIXLN) GO TO 2100
      LASTPX = NPIXLN
      MINNC = NPIXLN - S1 + 1
 2100 CONTINUE
C
C     NOW USE NUMBER OF INPUT COLUMNS AVAILABLE TO SET CELL SEGMENT
C     INDICATORS AND TO FIND ROW AND COLUMN INFORMATION
      MINS1 = 1000000
C
C     LOOP OVER CELL ROWS
      DO 2150 I=1,NHGL1
      IC1 = 100
      IC2 = 0
      RSTRT(I,NCGP) = 0
C
          DO 2130 J=1,NVGL1
C
C         CHECK WHETHER THE CELL REQUIRES INPUT OUTSIDE OF THE IMAGE
          IF (CELLCG(I,J).EQ.0) GO TO 2126
C
C         JUMP OUT IF THE CELL REQUIRES TOO MUCH CORE.
          IF (CELLCG(I,J).EQ.255) GO TO 2140
C
C         IF THE CELL USES INPUT DATA, CHECK FOR THE LAST PIXEL
          IF (CELLCG(I,J).EQ.222) GO TO 2125
          GO TO 2130
C
C         JUMP OUT IF THE CELL REQUIRES SAMPLES PAST LASTPX.
 2125     ROWNE = GRIDSL(2,I,J+1) - 1.0
          PIXNE = GRIDSL(1,I,J+1) + DELAY(ROWNE) + 3.0
          ROWSE = GRIDSL(2,I+1,J+1) + 2.0
          PIXSE = GRIDSL(1,I+1,J+1) + DELAY(ROWSE) + 3.0
          SE = MIN0 (MAX0(PIXNE,PIXSE), NPIXLN)
          IF (SE.GT.LASTPX) GO TO 2140
C
C         FIND FIRST AND LAST CELLS THAT USE INPUT DATA
          IC1 = MIN0 (J, IC1)
          IC2 = MAX0 (J, IC2)
C
C         ELSE SET COLUMN GROUP INDICATOR TO GROUP NUMBER.
          CELLCG(I,J) = NCGP
          GO TO 2130
C
 2126     CELLCG(I,J) = NCGP + 100
 2130     CONTINUE
C
C         FIND FIRST INPUT ROW FOR THIS ROW OF CELLS
 2140     IF (IC2.EQ.0) GO TO 2150
          ROWNW = GRIDSL(2,I,IC1) - 1.0
          ROWNE = GRIDSL(2,I,IC2+1) - 1.0
          RN = MAX0 (MIN0(ROWNW,ROWNE), 1)
          RSTRT(I,NCGP) = RN
C
C         FIND STARTING INPUT COLUMN FOR NEXT SEGMENT
          IF (J.EQ.NVGL) GO TO 2150
          NEXTCL = IC2 + 1
          ROWNW = GRIDSL(2,I,NEXTCL) - 1.0
          PIXNW = GRIDSL(1,I,NEXTCL) + DELAY(ROWNW) - 2.0
          ROWSW = GRIDSL(2,I+1,NEXTCL) + 2.0
          PIXSW = GRIDSL(1,I+1,NEXTCL) + DELAY(ROWSW) - 2.0
          SW = MAX0 (MIN0(PIXNW,PIXSW), 1)
          MINS1 = MIN0 (SW, MINS1)
C
C         SET INDICATOR FOR STARTING CELL IN NEXT SEGMENT
          CELLCG(I,NEXTCL) = 222
 2150 CONTINUE
C
C     SET PARAMETERS OF INPUT DATA SEGMENT
      MAXNR = MIN0 (CORE/MINNC, NLINPX)
      NBSAM(NCGP) = MINNC
      NBLIN(NCGP) = MAXNR
      NSTRT(NCGP) = S1
      S1 = MINS1
C
C     IF ALL CELLS NOT DONE, START A NEW SEGMENT
      IF (LASTPX.LT.NPIXLN) GO TO 1800
C
      WRITE (6,100)
      DO 2160 I=1,NHGL1
 2160 WRITE (6,101) (CELLCG(I,J), J=1,NVGL1)
      RETURN
C
  100 FORMAT (///5X,'COLUMN GROUP NUMBERS FOR THE GRID CELLS'/5X,
     .39('*')/)
  101 FORMAT (29(1X,I2))
      END
C
      SUBROUTINE RECTFI (INPIX, PIXEL, PIXOUT)
C
C      *CALL RESAMPLING ROUTINE TO COMPUTE INPUT PIXEL LOCATIONS AND
C       PERFORM SPECIFIED INTERPOLATION METHOD.
C       *READ LINE SEGMENTS AND WRITE ASSEMBLED RECORDS.
C
C     GRID IS COORDINATES (UTM & PIXEL) OF HORIZONTAL & VERTICAL GRID
C     INTERSECTIONS
C     UTMPIX IS THE FOLLOWING UTM TO PIXEL TRANSFORMATION PARAMETERS:
C     1. PIXEL COORDINATE OF FIRST SAMPLE IN SEGMENT
C     2. LINE COORDINATE OF FIRST SAMPLE IN SEGMENT
C     3. SLOPE OF INPUT PIXEL VS. OUTPUT PIXEL EQUATION
C     4. SLOPE OF INPUT  LINE VS. OUTPUT PIXEL EQUATION
C     5. NUMBER OF OUTPUT PIXELS IN THE CELL
C     6. OUTPUT PIXEL NUMBER.
C
      LOGICAL*1 INPIX(1), PIXEL(1), PIXOUT(NB,NSPL), CELLCG(29,29)
      DIMENSION GRIDEN(2,30), GRIDSL(2,30,30), NUMPXC(30), CPXVAL(30),
     .NUMLNC(30)
      INTEGER SAMPSG(29,10), RSTRT, PSTRT, SAMP0
      INTEGER*2 MNDEX(1000)
      REAL NORTH
      COMMON /PIXDAT/ NLINPX, NPIXLN, NHGL, NVGL, SRATE, NLPI, NSPL,
     .                KDUM, INTERP, NB
      COMMON /  GRID/ GRIDEN, GRIDSL, NUMPXC, CPXVAL, NUMLNC, CELLCG
      COMMON /SEGMNT/ NCGP, NBSAM(10), NBLIN(10), PSTRT(10),
     .                RSTRT(29,10)
      COMMON / PIXIN/ MAXSP, MAXLN, MNDEX, SAMP0, LINE0
      COMMON / TRANS/ UTMPIX(6)
      EQUIVALENCE (NPIX, UTMPIX(5)), (ISAMP, UTMPIX(6))
C
      NHGL1 = NHGL - 1
      NVGL1 = NVGL - 1
      IF (GRIDEN(2,2).GT.GRIDEN(2,1)) SRATEN =  SRATE
      IF (GRIDEN(2,2).LT.GRIDEN(2,1)) SRATEN = -SRATE
C
C     LOOP OVER COLUMN GROUPS
      DO 1020 ICGP=1,NCGP
      MAXSP = NBSAM(ICGP)
      MAXLN = NBLIN(ICGP)
      SAMP0 = PSTRT(ICGP)
      LUNIT= ICGP + 20
      REWIND LUNIT
      CALL LOADLN (.TRUE., 1, MAXLN, INPIX, PIXEL)
C
C     LOOP OVER OUTPUT IMAGE CELLS
      DO 1010 IGY1=1,NHGL1
      IGY2 = IGY1 + 1
      LSTRT = NUMLNC(IGY1)
      LSTOP = NUMLNC(IGY2) - 1
C
C       LOAD LINES OF CCT NECESSARY FOR RESAMPLED LINE
        LMIN = RSTRT(IGY1,ICGP)
        IF (LMIN.EQ.0) GO TO 150
        LMAX = MIN0 (LMIN+MAXLN-1, NLINPX)
        CALL LOADLN (.FALSE., LMIN, LMAX, INPIX, PIXEL)
C
C     LOOP OVER OUTPUT IMAGE LINES IN THE CELL
  150 DO 1000 ILINE=LSTRT,LSTOP
      NORTH = GRIDEN(2,1) + (ILINE-1)*SRATEN
C
C     LOOP OVER GRID CELLS COVERING MAP AREA
      ISAMP=0
      DO 300 IGX1=1,NVGL1
      IGX2 = IGX1 + 1
C
C           FIND NUMBER OF OUTPUT PIXELS IN THE CELL
            NPIX = NUMPXC(IGX2) - NUMPXC(IGX1)
C
C           CHECK FOR CELLS NOT REQUIRING INPUT DATA
            IF (CELLCG(IGY1,IGX1).EQ.ICGP+100) GO TO 250
            IF (CELLCG(IGY1,IGX1).EQ.ICGP) GO TO 210
            GO TO 300
C
C     LOAD UTM COORDINATES OF GRID CELL CORNERS
  210 AX = GRIDEN(1,IGX1)
      BX = GRIDEN(1,IGX2)
      AY = GRIDEN(2,IGY1)
      CY = GRIDEN(2,IGY2)
C
C     LOAD PIXEL COORDINATES OF GRID CELL CORNERS
      AS = GRIDSL(1,IGY1,IGX1)
      AL = GRIDSL(2,IGY1,IGX1)
      BS = GRIDSL(1,IGY1,IGX2)
      BL = GRIDSL(2,IGY1,IGX2)
      CS = GRIDSL(1,IGY2,IGX1)
      CL = GRIDSL(2,IGY2,IGX1)
      DS = GRIDSL(1,IGY2,IGX2)
      DL = GRIDSL(2,IGY2,IGX2)
C
C     COMPUTE FRACTIONAL CELL DISTANCE ALONG NORTHING AXIS
      ACRATO = (AY-NORTH) / (AY-CY)
C
C     COMPUTE INTERPOLATED SAMPLE AND LINE VALUES AT CELL EDGES
      QS = AS + ACRATO*(CS-AS)
      QL = AL + ACRATO*(CL-AL)
      RS = BS + ACRATO*(DS-BS)
      RL = BL + ACRATO*(DL-BL)
C
C     COMPUTE SLOPE OF INPUT SAMPLES AND LINES WITHIN THE CELL
      DELTAE = BX - AX
      UTMPIX(3) = (RS-QS) / DELTAE
      UTMPIX(4) = (RL-QL) / DELTAE
C
C           FIND INPUT COORDINATES OF FIRST POINT IN THIS GRID CELL
            UX1 = CPXVAL(IGX1) - AX
            UTMPIX(1) = QS + UX1*UTMPIX(3)
            UTMPIX(2) = QL + UX1*UTMPIX(4)
C
C           INTERPOLATE OVER INPUT DATA
            CALL RESAMP (PIXEL, PIXOUT)
            GO TO 300
C
C           FILL IN THE CELL NOT REQUIRING INPUT DATA
  250       DO 255 I=1,NPIX
            ISAMP = ISAMP + 1
            DO 255 J=1,NB
  255       PIXOUT(J,ISAMP) = KDUM
  300       CONTINUE
C
C         WRITE THE LINE SEGMENT OUT
          IF (NCGP.EQ.1) GO TO 500
          IF (ISAMP.EQ.0) GO TO 1010
          CALL WRITAR (LUNIT, PIXOUT, NB*ISAMP)
          GO TO 1000
C
  500     CONTINUE
          WRITE (11) PIXOUT
 1000     CONTINUE
C
 1010   SAMPSG(IGY1,ICGP) = ISAMP
C
      REWIND LUNIT
      WRITE (6,1700) LSTOP, ICGP
 1020 CONTINUE
C
C     ASSEMBLE AND WRITE OUTPUT IMAGE
      IF (NCGP.EQ.1) RETURN
      DO 1050 IGY1=1,NHGL1
      LSTRT = NUMLNC(IGY1)
      LSTOP = NUMLNC(IGY1+1) - 1
C
C     ASSEMBLE AND WRITE OUT THE LINE SEGMENTS
      DO 1040 ILINE=LSTRT,LSTOP
      NS = 1
      DO 1030 ICGP=1,NCGP
      NSAMP = SAMPSG(IGY1,ICGP)
      IF (NSAMP.EQ.0) GO TO 1030
      LUNIT= ICGP + 20
      CALL READAR (LUNIT, PIXOUT(1,NS), NB*NSAMP)
      NS = NS + NSAMP
 1030 CONTINUE
 1040 WRITE (11) PIXOUT
 1050 CONTINUE
      RETURN
C
1700  FORMAT(' FINISHED PROCESSING'I5,' RECORDS IN COLUMN GROUP'I3)
      END
C
      SUBROUTINE LOADLN (FIRST, LMIN, LMAX, INPIX, PIXEL)
C
C     THIS SUBROUTINE LOADS ARRAY 'PIXEL' WITH INPUT IMAGE LINES
C     LMIN THROUGH LMAX. IT FIRST CHECKS WHICH LINES ARE ALREADY
C     LOADED (L1 THROUGH L2) TO DETERMINE WHICH LINES CAN BE LEFT
C     IN MEMORY AND WHICH MUST BE READ.
C     LNDEX(I) IS THE LINE NUMBER OF DATA STORED IN PIXEL(*,I)
C     MNDEX(L) IS STORAGE LOCATION OF LINE NUMBER 'L'
C     IPOSN IS THE LINE NUMBER TAPE IS POSITIONED AT
C
      LOGICAL*1 INPIX(NB,NPIXLN), PIXEL(NB,NBSAM,MAXLN)
      INTEGER*2 LNDEX(1000),MNDEX(1000), SAMP0*4
      LOGICAL EOF, FIRST
      COMMON /PIXDAT/ NLINPX, NPIXLN, NHGL, NVGL, SRATE, NLPI, NSPL,
     .                KDUM, INTERP, NB
      COMMON / PIXIN/ NBSAM, MAXLN, MNDEX, SAMP0, LINE0
C
C     INITIALIZE BY CALLING LAST LINE IN CORE LINE '0'
      IF (.NOT.FIRST) GO TO 200
      NBSAM4 = NB*NBSAM
      REWIND 10
      EOF = .FALSE.
      IPOSN = 1
      L2 = MAXLN
      L1 = 1
      LNDEX(L1)=1-MAXLN
      LNDEX(L2)=0
      GO TO 1200
C
C     JFIRST, JLAST - FIRST AND LAST LINES ALREADY LOADED
  200 IF (EOF) RETURN
      JFIRST=LNDEX(L1)
      JLAST=LNDEX(L2)
      IF(LMIN.GE.JFIRST) GO TO 211
      WRITE(6,190) LMIN,LMAX,JFIRST,JLAST
      STOP 41
C
C     NLINRD - NUMBER OF LINES REQUIRED TO FILL IN FROM JLAST TO LMAX
C     JNEW   - NEW FIRST LINE AFTER LOADING NLINRD LINES FROM JFIRST
211   CONTINUE
      NLINRD=LMAX-JLAST
      JNEW=JFIRST+NLINRD
      IF(NLINRD.LE.0) GO TO 1200
      IF(NLINRD.GE.MAXLN) GO TO 300
C
C     FIND INDEX OF LINE JNEW (I)
C     LMINX1 - INDEX OF LAST LINE TO BE READ IN (I-1)
      DO 210 I=1,MAXLN
      IF(LNDEX(I).NE.JNEW) GO TO 210
      LMINX1 = I - 1
      IF(LMINX1.LE.0)LMINX1=MAXLN
      GO TO 220
  210 CONTINUE
  220 CONTINUE
      IF(L1.GT.LMINX1) GO TO 240
C
C *** CASE 1 - LOAD DATA FROM L1 TO LMINX1
      DO 230 I=L1,LMINX1
      READ (10,END=1300,ERR=1400) INPIX
      CALL MVL (INPIX(1,SAMP0), PIXEL(1,1,I), NBSAM4)
      LNDEX(I)=IPOSN
      IPOSN = IPOSN + 1
  230 CONTINUE
      L1 = LMINX1 + 1
      L2 = LMINX1
      IF(L1.GT.MAXLN) L1 =1
      GO TO 1000
C
C *** CASE 2 - LOAD DATA FROM L1 THROUGH MAXLN AND 1 THROUGH LMINX1
  240 DO 250 I = L1,MAXLN
      READ (10,END=1300,ERR=1400) INPIX
      CALL MVL (INPIX(1,SAMP0), PIXEL(1,1,I), NBSAM4)
      LNDEX(I) = IPOSN
      IPOSN = IPOSN + 1
  250 CONTINUE
      DO 260 I = 1,LMINX1
      READ (10,END=1300,ERR=1400) INPIX
      CALL MVL (INPIX(1,SAMP0), PIXEL(1,1,I), NBSAM4)
      LNDEX(I) = IPOSN
      IPOSN = IPOSN + 1
  260 CONTINUE
      L1 = LMINX1 + 1
      L2 = LMINX1
      GO TO 1000
C
C *** CASE 3 - NO OVERLAP OF OLD AND NEW DATA
C     POSITION TAPE
  300 IF(IPOSN.EQ.JNEW)GO TO 320
      READ (10,END=1300,ERR=1400)
      IPOSN =IPOSN + 1
      GO TO 300
  320 DO 330 I=1,MAXLN
      READ (10,END=1300,ERR=1400) INPIX
      CALL MVL (INPIX(1,SAMP0), PIXEL(1,1,I), NBSAM4)
      LNDEX(I) = IPOSN
      IPOSN = IPOSN + 1
  330 CONTINUE
      L1 = 1
      L2  = MAXLN
C
C     LOAD STORAGE LOCATIONS INTO ELEMENT OF MNDEX EQUAL TO LINE NUMBER
C     LINE0 - FIRST IMAGE LINE HELD IN CORE
C     MX    - STORAGE LINE NUMBER AT IMAGE LINE LNDEX(I)
 1000 LINE0 = LNDEX(L1)
      DO 1100 I = 1,MAXLN
      MX = LNDEX(I) - LINE0 + 1
      MNDEX(MX) = I
 1100 CONTINUE
 1200 CONTINUE
      RETURN
C
 1300 EOF = .TRUE.
      WRITE(6,1301) IPOSN
      STOP 41
 1400 WRITE (6,1401) IPOSN
      STOP 43
C
190   FORMAT(1H0,'ERROR: BACKWARD READ REQUESTED',4I10)
1301  FORMAT(1H0,'EOF AT LINE ',I5)
1401  FORMAT(1H0,'READ ERROR AT LINE',I5)
      END
C
      SUBROUTINE RESAMP (LPIXEL, LPXOUT)
C
C     INTERP = 1, 2, 3  GIVES NEAREST, BILINEAR, BICUBIC INTERPOLATION
C
      LOGICAL*1 LPIXEL(NB,NBSAM,MAXLN), LPXOUT(NB,1)
      DIMENSION RK(4,4)
      INTEGER SAMP0, SLR
      INTEGER*2 MNDEX(1000)
C
C     UTMPIX IS THE FOLLOWING UTM TO PIXEL TRANSFORMATION PARAMETERS:
C     1. PIXEL COORDINATE OF FIRST SAMPLE IN SEGMENT
C     2. LINE COORDINATE OF FIRST SAMPLE IN SEGMENT
C     3. SLOPE OF INPUT PIXEL VS. OUTPUT PIXEL EQUATION
C     4. SLOPE OF INPUT  LINE VS. OUTPUT PIXEL EQUATION
C     5. NUMBER OF OUTPUT PIXELS IN THE CELL
C     6. OUTPUT PIXEL NUMBER.
      COMMON /PIXDAT/ NLINPX, NPIXLN, NHGL, NVGL, SRATE, NLPI, NSPL,
     .                KDUM, INTERP, NB
      COMMON / PIXIN/ NBSAM, MAXLN, MNDEX, SAMP0, LINE0
      COMMON / TRANS/ UTMPIX(6)
      EQUIVALENCE (NPIX, UTMPIX(5)), (ISAMP, UTMPIX(6))
C
      LBASE=LINE0-1
      PSI = UTMPIX(1) - SAMP0 + 1
      PLI = UTMPIX(2) - LINE0 + 1
      SPACE1 = UTMPIX(3)*SRATE
      SPACE2 = UTMPIX(4)*SRATE
C
C         LOOP OVER OUTPUT PIXELS
          DO 1050 I=1,NPIX
          ISAMP = ISAMP + 1
          DE=I-1
C
C         COMPUTE INPUT LINE AND SAMPLE NUMBERS
          PS=PSI+DE*SPACE1
          PL=PLI+DE*SPACE2
          GO TO (1000, 2000, 3000), INTERP
C
C     * * * * * NEAREST NEIGHBOR INTERPOLATION * * * * *
C
 1000     IPL = PL + 0.5
          IF (IPL.LT.1.OR.IPL.GT.MAXLN) GO TO 1030
C
C           ADD EARTH ROTATION AND SENSOR DELAY OFFSETS
            IPS = PS + DELAY(IPL+LBASE) + 0.5
            IF (IPS.LT.1.OR.IPS.GT.NBSAM) GO TO 1030
            SLR = MNDEX(IPL)
C
          DO 1021 IBAND=1,NB
 1021     LPXOUT(IBAND,ISAMP) = LPIXEL(IBAND,IPS,SLR)
          GO TO 1050
C
C     * * * * * BILINEAR INTERPOLATION * * * * *
C
 2000     IPL = PL
          IF (IPL.LT.1.OR.IPL.GT.MAXLN-1) GO TO 1030
          DP = PL - IPL
            DO 2010 LINE=1,2
            PSCCT = PS + DELAY(IPL+LBASE)
            IPS = PSCCT
            IF (IPS.LT.1.OR.IPS.GT.NBSAM-1) GO TO 1030
            D = PSCCT - IPS
            SLR = MNDEX(IPL)
              DO 2005 IBAND=1,NB
              K1 = LPIXEL(IBAND,IPS  ,SLR)
              K2 = LPIXEL(IBAND,IPS+1,SLR)
 2005         RK(LINE,IBAND) = K1 + D*(K2-K1)
            IPL = IPL + 1
 2010       CONTINUE
          DO 2021 IBAND=1,NB
          IX1 = RK(1,IBAND) + DP*(RK(2,IBAND)-RK(1,IBAND)) + 0.5
 2021     LPXOUT(IBAND,ISAMP) = IX1
          GO TO 1050
C
C     * * * * * BICUBIC INTERPOLATION * * * * *
C
 3000     IPL = PL
          IF (IPL.LT.2.OR.IPL.GT.MAXLN-2) GO TO 1030
          DP = PL - IPL
          IPL = IPL - 1
C
C           COMPUTE AN INTERPOLATED SAMPLE IN EACH OF 4 LINES
            DO 3010 LINE=1,4
            PSCCT = PS + DELAY(IPL+LBASE)
            IPS = PSCCT
            IF (IPS.LT.2.OR.IPS.GT.NBSAM-2) GO TO 1030
            D = PSCCT - IPS
            SLR = MNDEX(IPL)
              DO 3005 IBAND=1,NB
              K1 = LPIXEL(IBAND,IPS-1,SLR)
              K2 = LPIXEL(IBAND,IPS  ,SLR)
              K3 = LPIXEL(IBAND,IPS+1,SLR)
              K4 = LPIXEL(IBAND,IPS+2,SLR)
              K21 = K2 - K1
              K43 = K4 - K3
 3005         RK(LINE,IBAND) = D* (D* (D* (K43+K21) - (K43+2*K21))
     .                         + (K3-K1)) + K2
            IPL = IPL + 1
 3010       CONTINUE
C
C         INTERPOLATE OVER 4 LINES TO GET FINAL OUTPUT SAMPLE
          DO 3021 IBAND=1,NB
          R21 = RK(2,IBAND) - RK(1,IBAND)
          R43 = RK(4,IBAND) - RK(3,IBAND)
          IX1 = DP* (DP* (DP* (R43+R21) - (R43+2.0*R21)) + (RK(3,IBAND)-
     .          RK(1,IBAND))) + RK(2,IBAND) + 0.5
          IF (IX1.LT.0) IX1=0
          IF (IX1.GT.255) IX1=255
 3021     LPXOUT(IBAND,ISAMP) = IX1
          GO TO 1050
C
 1030 DO 1035 IBAND=1,NB
 1035 LPXOUT(IBAND,ISAMP) = KDUM
C
 1050     CONTINUE
      RETURN
      END
C
      SUBROUTINE EVPOLY (IFUN, X, Y, ANS)
C
C     EVALUATE POLYNOMIAL FIT FUNCTIONS
C
      DOUBLE PRECISION COEF(21,4), ANSD, XD, YD, XXX, YYY
      INTEGER NTERM(5)/3,6,10,15,21/
      INTEGER XP(21)/0,1,0,2,1,0,3,2,1,0,4,3,2,1,0,5,4,3,2,1,0/
      INTEGER YP(21)/0,0,1,0,1,2,0,1,2,3,0,1,2,3,4,0,1,2,3,4,5/
      COMMON /LSQCFC/ COEF, LSQDEG, IER1(2), TOL
      COMMON / MEANS/ GCPM(4)
C
C     COEF(21,I)- COEFFICIENTS FOR : 1- PIXEL
C                                    2- LINE
C                                    3- EASTING
C                                    4- NORTHING
      ISTOP = NTERM(LSQDEG)
      ANSD = 0.0
      IF (IFUN.EQ.3.OR.IFUN.EQ.4) GO TO 1
      XD = X - GCPM(3)
      YD = Y - GCPM(4)
      GO TO 5
    1 XD = X - GCPM(1)
      YD = Y - GCPM(2)
    5 CONTINUE
C
      DO 10 I=1,ISTOP
      XXX = 1.0
      IF (XP(I).NE.0) XXX = XD**XP(I)
      YYY = 1.0
      IF (YP(I).NE.0) YYY = YD**YP(I)
   10 ANSD = ANSD + COEF(I,IFUN)*XXX*YYY
C
      ANS = ANSD + GCPM(IFUN)
      RETURN
      END
C
      FUNCTION DELAY (LINE)
C
C     COMPUTE ROTATIONAL AND SENSOR DELAY RELATIVE TO FIRST SWATH
C
C     DEGCEN - LATITUDE AT THE CENTER OF THE LANDSAT SCENE
C     PIXDLY - EQUATORIAL EARTH ROTATION PER SWATH IN PIXELS
C     DEGPLN - CHANGE IN LATITUDE PER SCAN LINE
C     RADDEG - RADIANS PER DEGREE
C     SDELAY - SENSOR SAMPLING INTERVAL BETWEEN LINES (2) / TOTAL NUMBER
C              OF DETECTORS (25)
C     LINESW - LINE NUMBER IN THE SWATH (1 - 6)
C
      DIMENSION SDELAY(6)
      INTEGER CCTLIN, SWATH
      LOGICAL CCT
      COMMON /LANDST/ CCT, LLC, DEGCEN, SAMPOF, LINOFF, AMPL, PHASE
      DATA PIXDLY, DEGPLN, RADDEG, SDELAY /0.6, 0.00071086, 0.01745329,
     .0.0, 0.08, 0.16, 0.24, 0.32, 0.40/
C
      IF (CCT) GO TO 10
      DELAY=0.0
      RETURN
C
   10 CONTINUE
      CCTLIN = LINE + LINOFF
      DEGLAT = DEGCEN + (1170.5-CCTLIN) * DEGPLN
      RDELAY = PIXDLY * COS (DEGLAT*RADDEG)
      SWATH = (CCTLIN-1)/6 + 1
      LINESW = MOD(CCTLIN-1,6) + 1
      DELAY = RDELAY*(SWATH-1) - SDELAY(LINESW)
      RETURN
      END
C
      REAL FUNCTION MVPOFF (PS)
C
C     COMPUTE MIRROR VELOCITY PROFILE OFFSET
C
C     LLC - NUMBER OF PIXELS IN THE RAW SCAN LINE
C     SAMPOF - NUMBER OF PIXELS SKIPPED IN THE LANDSAT SCENE
C     AMPL, PHASE - AMPLITUDE, PHASE OF MIRROR VELOCITY PROFILE CURVE
C
      LOGICAL CCT
      COMMON /LANDST/ CCT, LLC, DEGCEN, SAMPOF, LINOFF, AMPL, PHASE
C
      IF (CCT) GO TO 12
      MVPOFF = 0.0
      RETURN
C
   12 CONTINUE
      PS1 = PS + SAMPOF
      MVPOFF = AMPL * SIN (6.2831853 * (PS1+PHASE-1.0) / (LLC-1))
      RETURN
      END
C
      REAL FUNCTION MVPINV (PS)
C
C     COMPUTE NEGATIVE MIRROR VELOCITY PROFILE CORRECTION
C
C     LLC - NUMBER OF PIXELS IN THE RAW SCAN LINE
C     SAMPOF - NUMBER OF PIXELS SKIPPED IN THE LANDSAT SCENE
C     AMPL, PHASE - AMPLITUDE, PHASE OF MIRROR VELOCITY PROFILE CURVE
C
      LOGICAL CCT
      COMMON /LANDST/ CCT, LLC, DEGCEN, SAMPOF, LINOFF, AMPL, PHASE
C
      IF (CCT) GO TO 10
      MVPINV = 0.0
      RETURN
C
   10 CONTINUE
      PS1 = PS + SAMPOF
      MVPINV = AMPL * SIN (6.2831853 * (PS1+PHASE-1.0) / (LLC-1))
      RETURN
      END
C
      FUNCTION ERCURV (PS)
C
C     COMPUTE EARTH CURVATURE CORRECTION
C
C     RE, RSAT - EARTH RADIUS, SATELLITE ORBIT RADIUS
C     TFOV - TOTAL FIELD OF VIEW OF THE SCANNER (11.56 DEGREES)
C     LLC  - NUMBER OF PIXELS IN THE RAW SCAN LINE
C     SAMPOF - NUMBER OF PIXELS SKIPPED IN THE LANDSAT SCENE
C
      LOGICAL CCT
      COMMON /LANDST/ CCT, LLC, DEGCEN, SAMPOF, LINOFF, AMPL, PHASE
      DATA RE, RSAT, TFOV /6367.4, 7285.6, 0.20176/
C
      IF (CCT) GO TO 10
      ERCURV = 0.0
      RETURN
C
C     COMPUTE SCANNER ANGLE AT PIXEL NO. PS AND ANGLE SUBTENDED AT
C     THE CENTER OF THE EARTH
   10 CONTINUE
      TFOV2 = TFOV/2.0
      PS1 = PS + SAMPOF
      ANGSCN = (PS1-1.0)*TFOV/(LLC-1) - TFOV2
      ANGERT = ARSIN(SIN(ANGSCN)*RSAT/RE) - ANGSCN
      TOTERT = ARSIN (SIN(TFOV2)*RSAT/RE) - TFOV2
C
C     FIND SCANNER ANGLE BASED ON FRACTION OF TOTAL ARC ON THE EARTH'S
C     SURFACE AND CONVERT TO PIXEL NUMBER
      ANGSC1 = TFOV2 * ANGERT / TOTERT
      PS2 = 1.0 + (ANGSC1+TFOV2) * (LLC-1) / TFOV
      ERCURV = PS1 - PS2
      RETURN
      END
C
      FUNCTION ERCURI (PS)
C
C     COMPUTE EARTH CURVATURE DE-CORRECTION
C
C     RE, RSAT - EARTH RADIUS, SATELLITE ORBIT RADIUS
C     TFOV - TOTAL FIELD OF VIEW OF THE SCANNER (11.56 DEGREES)
C     LLC  - NUMBER OF PIXELS IN THE RAW SCAN LINE
C     SAMPOF - NUMBER OF PIXELS SKIPPED IN THE LANDSAT SCENE
C
      LOGICAL CCT
      COMMON /LANDST/ CCT, LLC, DEGCEN, SAMPOF, LINOFF, AMPL, PHASE
      DATA RE, RSAT, TFOV /6367.4, 7285.6, 0.20176/
C
      IF (CCT) GO TO 10
      ERCURI = 0.0
      RETURN
C
C     EARTH ANGLE IS PROPORTIONAL TO CORRECTED PIXEL NUMBER.
   10 CONTINUE
      TFOV2 = TFOV/2.0
      CENTER = (LLC+1)/2.0
      PS1 = PS + SAMPOF
      TOTERT = ARSIN (SIN(TFOV2)*RSAT/RE) - TFOV2
      ANGERT = -TOTERT * (CENTER-PS1) / (CENTER-1.0)
C
C     COMPUTE ORIGINAL SCAN ANGLE BASED ON EARTH ANGLE.
C     C IS LINE FROM SATELLITE TO PIXEL LOCATION ON THE EARTH'S SURFACE.
      C = (RSAT**2 + RE**2 - 2.0*RSAT*RE*COS(ANGERT)) ** 0.5
      ANGSCN = ARSIN(RE/C*SIN(ANGERT))
      PS2 = 1.0 + (ANGSCN+TFOV2) * (LLC-1) / TFOV
      ERCURI = PS2 - PS1
      RETURN
      END
C
      SUBROUTINE READAR (NTAPE1, W, NSAMP)
C
C  READ NSAMP BYTES INTO ARRAY W FROM LOGICAL UNIT NTAPE1
C
      LOGICAL*1 W(NSAMP)
C
      READ (NTAPE1) W
      RETURN
C
C     ******************************************************************
C
      ENTRY WRITAR (NTAPE1, W, NSAMP)
C
C  WRITE NSAMP BYTES FROM ARRAY W ONTO LOGICAL UNIT NTAPE1
C
      WRITE (NTAPE1) W
      RETURN
      END
